clear;clc
close all

%% init arguments
Fs = 1e16;  % sampling frequence /Hz
T = 1/Fs;   % sampling time /s
Number0 = 1000000; % original signal length
Number = 1000000;  % data length
t0 = (0:Number0-1)*T;   % original signal time sequence
t = (0:Number-1)*T; % data time sequence
y = zeros(1,Number);    % init 

light_labmda = 5.32*10e-8;    
sound_labmda =  1.5*10e-4;      
n0 = 1.33;             
L = (sound_labmda^2*n0)/(2*pi*light_labmda);          
Deltan = 1.26*10e-12*101325;            
c = 3e+9;             
v = 1500;              
light_ommiga = c/light_labmda;         
sound_ommiga = v/sound_labmda;             
K = 2*pi/sound_labmda;        
Phi0 = (2*pi/light_labmda)/(n0*L);
Phi1 = (2*pi/light_labmda)/(Deltan*L);
Z = 0;
J0 = besselj(0,Phi0);  
J1 = besselj(1,Phi1);  
Et = J0 *exp((light_ommiga*t0 - Phi0)*1i) - J1 *(exp(((sound_ommiga + light_ommiga)*t0 - K*Z - Phi0)*1i) - exp(((light_ommiga - sound_ommiga)*t0 + K*Z - Phi0)*1i));
y(1:Number0) = Et;

%% plot

%plot(t,abs(Et))
%fourier_Et=fft(Et); 
%plot(f,abs(fourier_Et))

%% FFT
Y = fft(y);
P2 = abs(Y/Number0);
P1 = P2(1:Number/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(Number/2))/Number;   % why
plot(f,P1)
axis([8.9648e14 8.9800e14 0 0.6e-4])    % lock dispaly area
title('Single-Sided Amplitude Spectrum of S(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
